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The problem of formulating self-consistent local and global 
stability exponents is shown to require global separation of 
variables. Posing the separation of variable problem, we see 
that many such separations are possible, but only one is con- 
sistent with both Hamiltonian dynamics and the boundedness 
requirement for a Lyapunov transform: the determinant of 
the modal matrix must be constant. Such stability exponents 
are invariant to any linear transformation of variables, and 
both the local stability exponents and modal matrix appear 
to be point functions in the original space, and introduce a 
true coordinate frame. Methods are presented to perform this 
separation at equlibrium points, about periodic orbits, and 
along general trajectories. Results of numerical experiments 
are given. 
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I. INTRODUCTION 



A nonlinear dynamical system can be written as 



X = f(X). 



(1) 



We shall consider only autonomous dynamical systems 
in this paper. If the system itself is time dependent, 
then defining an additional state variable xn + ± = t and 
appending xn+i = 1 to (Q) will render the system au- 
tonomous. Writing x = X(t) — Xo(f) as the difference 
between a nominal solution Xo and a nearby trajectory, 
to first order x obeys 



df 



x = A(i)x. 



(2) 



As a linear system, the complete solution to the above is 
contained within the fundamental matrix which obeys 



$ = A(t)$, $(*o) = I- 



(3) 



At least for finite final times, the fundamental matrix can 
be constructed by numerical integration. The general 
solution to (0) is then written as 



x(t) =$(t,to)x(t ), 



(4) 



where the two time indices on $ indicate the end and 
start times, respectively. 

The stability and control of chaotic systems are usu- 
ally discussed in terms of Lyapunov exponents. As is 



well known, Lorenz ||, Greene and Kim [Q, Lyapunov 
exponents appear in the singular value / singular vector 
decomposition of the fundamental matrix 



*(i/,*o) = Uexp(A(tf-to))V 7 



(5) 



where U and V are real orthonormal matrices, and where 
the elements of the real diagonal matrix A are the Lya- 
punov exponents. As tf — * oo, the orthonormal matrix 
V approaches a constant, and a local coordinate frame is 
introduced near the trajectory. 

But this decomposition does not match the usual de- 
compositions for constant coefficient and periodic coeffi- 
cient systems. For constant coefficient systems, the usual 
decomposition is 



$(t f) to)=Eexp(A(t f -to))E- 



(6) 



where E is the matrix of eigenvectors and A is the Jordan 
normal form containing the eigenvalues of the constant 
matrix A. This form is only possible because A is con- 
stant. In the periodic coefficient case, Floquet showed 
that the fundamental matrix decomposed as 



$(t, t ) = E{t) exp(A(t - to^E-^to). 



(7) 



Here, the matrix E(t) is time dependent but periodic, 
while the Jordan form A contains the Poincare exponents. 
The most obvious difference between (^), (|^) and (||) is 
that the former two permit imaginary parts in the stabil- 
ity exponents, while Lyapunov exponents are purely real 
quantities. 

There are more subtle differences as well between Lya- 
punov exponents and the other stability exponents. The 
Floquet decomposition, when applied to an equilibrium 
point immediately reduces to the usual constant coeffi- 
cient solution. The decomposition does not reduce 
to either of the other two cases. Also, the state vector 
X usually consists of disparate physical quantities, often 
not measured in the same physical units. Assume two 
realizations of the state vector are X and Y, related by 
the coordinate transformation law 



X{Y). 



(8) 



Then two representations of the local state vector are 
related by 



x(i) = M(t)y(t), 



(9) 



where M — dX /dY is the Jacobian matrix of the trans- 
formation of coordinates. [The Jacobian is actually a 



1 



function of position X, but we use the reference trajec- 
tory X (i) to write it as a function of time. This is, 
however, only a convenient shorthand. We are consid- 
ering autonomous coordinate transformations, and M is 
not an explicit function of time.] Using the above in 
(Q), and remembering that the fundamental matrix is 
the derivative of the final state with respect to the initial 
state, $ y (t, to) = dy(t)/dy(t ), the fundamental matrix 
expressed in the variables y is 



%(t f ,t )=M- 1 (t f )^ x {t f ,t )M(t ). 



(10) 



Inspection of (|) and (0) shows that these forms are 
compatible with this transformation. In particular, the 
matrix A of stability exponents is invariant to any au- 
tonomous change in coordinate systems. This is true in 
the constant coefficient case since the reference trajec- 
tory is a point, and so M(tf) — M(to). This makes 
similar to Q x . In the periodic case, the stability expo- 
nents are evaluated from the <!> matrix at one full period, 
and so again the Jacobian matrices are the same, and 
the $ matrices are similar. Stability exponents in these 
cases are fully invariant to any autonomous coordinate 
transformation. 

On the other hand, the classical Lyapunov exponents 
(|J) are only invariant under rotations of the original co- 
ordinate system. They are not even invariant under 
changes in physical units. This is a nontrivial matter 
when one is not working in an abstract "metric space", 
but instead has to contend with physical units for a phys- 
ical problem. 

Also, the "local" Lyapunov exponents deserve mention 
here, e.g. H, ||. They are the eigenvalues of A T A at a 
point, and are an attempt to define a quantity which 
would characterize the instantaneous growth and decay 
rates of small displacements. However, since they ignore 
terms relating to the rate of change of the eigenvectors of 
A T A, they do not integrate with time to give the global 
Lyapunov exponents. 

In earlier works we have explored extensions of the con- 
stant coefficient and periodic coefficient decompositions 
to the general case of (§). In Wiesel M, we used the 
decomposition $(t,t ) = E(t) exp(A(*X* - io))25 _1 (t), 
which is just the instantaneous eigenvalue decomposition. 
But the stability exponents of this form will not be in- 
variant to most coordinate transformations either, and 
the assumption that E(t) = E(t ) seems unnatural. In 
|pl| the decomposition (0) was extended to more general 
systems, and winding numbers were successfully calcu- 
lated as the imaginary parts of the constant matrix A. 
But that paper only investigated deterministic regions 
of conservative Hamiltonian systems. In this paper we 
will explore the conditions under which the invariance 
of the stability exponents of the constant coefficient case 
can be extended, and the extent to which self-consistent 
local and global stability exponents can be defined. 



II. LYAPUNOV TRANSFORMS 

Consider a slight extension of the Floquet decomposi- 
tion (0) to 



$(t/,tb) = E(t) expeno))^- 1 ^). 



(11) 



The matrix Q will be presumed diagonal, since it must 
commute with its derivative for what is to follow. Sta- 
bility exponents over the time interval then are 



(12) 



t f - 1 

Inserting ( O ) into (||) and rearranging produces 

E = AE - EQ, (13) 

which assumes that Q and Q(t) commute. If the matrix 
determinant of E(t) is bounded from both above and 
below, 



\\E(t)\\<p, WE-WW^p, 



(14) 

for all t within to < t < tf then ( |ll| ) is a Lyapunov 
transformation, Rugh ||. That is, it is an example of 
the most general type of linear transformation of (|^) that 
will preserve stability information. We will refer to the 
LOi [and their limits as i/ -> oo] as extended stability 
exponents. They may have imaginary parts. 

Without the conditions on the determinant of E(t), the 
transformation (11) is very general. If we specify any di- 
agonal matrix function f2(t) and choose a random E(to), 
then ( |I3| ) will produce an E(tf) which will reproduce 
<&(tf,to) through (|ll| ). I n most cases the determinant 
of E will collapse or expand exponentially. This is not 
desirable, since this matrix can be interpreted as a local 
coordinate transformation. Let 

x = £(i)y. (15) 

Then calculating x, and using @ and (£f), wc find 



y = O y . 



(16) 



That is, the transformation ( |l5| ) will decouple the equa- 
tions of motion of the linear system. This decoupling is 
only meaningful, however, if E and its inverse are well 
behaved, so that the coordinate transformation ( |l5| ) is 
legitimate. 

It is interesting, however, to compare ( |l6| ) to the anal- 
ogous form (U). Suppose that we calculated the varia- 
tional equations (||) for a dynamical system, and found 
that they were of the form (|l^)? Then (|l6| ) would not 
be a local diagonalization along one particular trajectory, 
but would diagonalize all trajectories simultaneously. If 
the new form (|l^) is globally the variational equations of 
a transformed version of the original dynamical system, 
then the system that gives rise to ([l6]) must have the 
form 



2 



gi(Yi). 



(17) 



The important thing to notice here is that the above 
equations, unlike ([!]), are decoupled. We now reverse the 
entire chain of deduction that led to this point. In order 
to produce local stability exponents fl which integrate 
along any trajectory to give global stability exponents 
f2, the dynamical system must be decoupled. The exis- 
tence of a global Lyapunov transformation implies such 
a decoupling, and vice versa. 



III. SEPARATION NEAR AN EQUILIBRIUM 

In this section we will concentrate on the vicinity of 
equlibrium points. Eq ( |l5| ) then becomes the linear por- 
tion of a series expansion about an equlibrium 

1 dE ia 
Xi = E ia y a + — — y a yp 



2! dyp 



1 d 2 E ia 
3! dypdy. 



(18) 



[Greek indices indicate summations, while roman indices 
are not summed. Also, instead of noting that these 
quantities are evaluated at an equlibrium, we will find 
it more convenient to explicitly show functional depen- 
dence wherever a quantity is not evaluated at the equlib- 
rium.] We note that not all such series expansions will 
yield true coordinate frames. The condition that the new 
y be a coordinate frame is symmetry in all indices except 
the first: 



8E U 



dE, 



03 



d 2 E u 



d 2 E, 



i/3 



d 2 E r 



dyp dy a ' dypdy 7 dy a dy 7 dypdy^ 



(19) 



and so forth. [The familiar Lie bracket conditions apply 
to the covariant derivatives, not the contravariant quan- 
tities above.] Along with this expansion, we have the 
parallel expansion of the new equations of motion about 
the equlibrium 



9i = ^iUi + 



1 dfl, 



■Vi 



i d 2 n 



3! dy: 



2 yf + 



(20) 



Separation of variables mandates that each Qi be a func- 
tion only of the corresponding . As an immediate corol- 
lary, each local stability exponent £li is constant on sur- 
faces of constant j/j. 

We begin by noting in a general coordinate transfor- 
mation X = X(Y) that 



Y = g = i?- 1 (Y)f(X(Y)), 



(21) 



where E = <9X/<9Y is the Jacobian of the transformation. 
At the equlibrium point, then, gi = 0. Proceeding to the 
first order, we have 



dy 3 



A 



Y,ij 



(22) 



E _ x df a | dE£ 
ia dy dyj 



— F-l j? it. 



i dEpy 



Eja fa 



E^Ax^pEpj - Ej-—^-E ^f a 



We have replaced the y partial of f with its equivalent in 
terms of x, recognized Ax = df/dx, and expanded the 
derivative of the inverse matrix. Then, evaluating at the 
equlibrium point, we must have 



A 



Y,ij 



fi = E~}A 



(23) 



To effect a separation of variables this must be diagonal, 
and we are immediately forced into using the eigenvalue / 
eigenvector decomposition of Ax as the first order values 
for Cli and E. 

At the second order, we calculate the second partial 
derivative of ( pl|) , and evaluate it at the equlibrium. The 
result is 



= El 



d 2 9i 

dy 3 dy k la [ dy k 

dAx,af3 



dE aj . _dE*k h 



dE, 



9Vj 



a 3 o 

dyk 3 



dx~, 



E-ykEpj 



(24) 



Again, for a separation of variables we must require that 
the above expression be zero, except for i = j = fc, 
when it must equal dtli/dyi. Note that ( p4| ) constitutes 
iV 3 linear equations in the N 3 unknowns dEij/dyk, but 
with the N additional unknown quantities d 2 gi/dy 2 — 
dQi/dyi. For the moment we delay specifying additional 
information to choose the quantities dfli/dyi. 

What we are attempting is similar to center manifold 
theory, and reduces to it if we drop all of the "diagonal" 
terms from (|2^), e.g. Arrowsmith and Place In center 
manifold theory only the surfaces containing the trajecto- 
ries are obtained. Our aim in this work goes beyond just 
obtaining the manifolds through the equlibrium. Also, 
this approach is similar to normal form theory. In fact, 
if we chose all the derivatives of fli to be zero, we would 
be attempting to map the entire dynamical system back 
onto the equlibrium point variables, and this would be 
normal form theory. This is also not our aim. 

Continuing to the third order we obtain 



dyjdykdyi 
dE^ 1 d 2 f a 



e: 



d 3 f a 



dyjdykdyi 
dEr 1 2 f a 



(25) 



dE- 1 d 2 f a 



dy.j dy k dyi dy k dy^dyi dyi dyjdy k 
d 2 E~ 1 df a | d 2 E- a 1 df a | d 2 E^ df a ^ 
dyjdyk dyi dy 3 dyi dy k dy k dyi dyj ' 



where 



3 



dy 3 



A,, B E, 



d 2 h 
dy dy k 



-E 7 kEfjj + Aifj- 



dy k ' 



(26) 



(27) 



d 2 A 



dyjdy k dyi dx^dxs 



i0 



EsiE 1 kEpj 



dA if3 ( dE, 



dx~, 



+ A U 



d 2 E, 



dyi 

0i 



E. 



dE, 



03 



dyi 



E. 



dE, 



03 



dyk 



E 



and 



dykdyi 



dykdyi 
_ 1 d 2 E i 



x dEjia_ ! 

* d Vj ° kl 



(28) 



(29) 



dE lS3 dE, 



0a 



dyi dy k 



E 



ip- 1 0a rp-l 

V dy k dyi °i 



E 



dyk dyi 



i0 



(30) 



Again, to force separation of variables, (|25|) must equal 
zero except for the N cases where i = j = k — I, in 
which case it equals d 2 (li/dy 2 . These are thus iV 4 linear 
equations in the N A + N unknowns d 2 Eij/ dykdyi and 
d^/dy 2 . 

Before proceeding to methods to specify the extra vari- 
ables at each order, we will first establish some properties 
of this transformation. First, gi = £'~^(Y)f Q (X(Y)) is 
presumably a continuous, differentiable function of Y, 
and therefore dgi/dyjdyk = dgi/dykdyj. Inserting ( gj ) 
into this expression and using the symmetry of deriva- 
tives of A, three terms immediately cancel, leaving 



E 



dyk 



(31) 



and a simple further simplification confirms the symme- 
try of dE a j j dyk- This is required if the new variables 
Y are to form a coordinate frame. At the third order, 
extensive numerical calculations have failed to produce a 
non-symmetric d 2 Eij/ dykdyi. 

To see why ( p4] ) and ( p5j ) leave the stability exponent 
unspecified, consider the simple dynamical system 



\J x 2 + y 2 



\Jx 2 +y- 



x-y, 



+ x-y. 



(32) 



This is the rectangular form of the polar variable differ- 
ential equations 



1 - r, 



1, 



(33) 



so this system obviously separates in polar coordinates 
Y\ = r, Y% = 9. But it also separates in any coordinate 
frame y which itself is a separable function of the polar 
coordinates 



yi = h 1 {r), y 2 = h 2 {6). 



(34) 



There is thus an infinite number of coordinate frames in 
which a separable system can be separated. The undeter- 
mined stability exponent derivatives in (|24|) and ( p5| ) are 
directly related to the derivatives of the arbitrary func- 
tions hi above. 

So, if the system separates in terms of the variables Y, 
then it also separates in any variables 3^ = D^C^i), where 
each yi is an arbitrary function of one Yi . The equations 
of motion in terms of these new variables take the form 



y - 1 s> 



(35) 



The local stability exponents for the new variables are 
&i,y = d^i /dyi, and direct calculation yields 



9., 



n 



i.Y 



d 2 Y ( dY 

^{y)g i {Y i {y i )) { ^(y) 



dy 



dy l 



(36) 



This result shows that these stability exponents are in- 
variant to constant coefficient linear transformations, 
where the second derivative is zero. If we study the same 
equlibrium beginning in two different coordinate frames, 
we will obtain the same first order stability exponents, 
and the eigenvectors will have the same direction, but al- 
most inevitably different magnitudes. The above result, 
however, shows that the higher order stability exponent 
derivatives will yield the same exponents when evaluated 
at the same point X in physical space. This result has 
also been confirmed numerically. It implies that Qi (X) is 
a true point function of the original space position vector 
X. [This is not true of f2, which for finite times will be 
a function of the arc studied.] This invariance class is 
also significantly stronger than the usual Lyapunov ex- 
ponents, which are only invariant under rigid rotations 
of the original coordinate frame. 

To ensure that the bound conditio ns (|l4| ) are met, and 
that the coordinate transformation (|15|) is nonsingular, 
we investigate the determinant of E(Ty Calculating the 
derivative of the determinant and using (^), one obtains 



JV 



-\E(t)\ =^|e 1 ...e i ...e w | (37) 

»=1 

N N 

= ^\e l ...AB i ...e rr \-^2(li\E{t)\. 



i=l 



The first term can then be reduced by a standard argu- 
ment M to produce 
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j t \E(t)\=(Tr(A(t))-Tr(n))\E(t)\, (38) 
where Tr() is the trace. This has solution 

\E(t)\ = \E(t )\ exp | jf (lV(A(t)) - Tt(fi(t))) dt 

(39) 

This form directly shows that \E(t) \ cannot become zero, 
so a lower bound exists over any finite time interval. Re- 
membering that the are local stability exponents, to 
obtain an upper bound in the long term it is most con- 
sistent to impose the instantaneous condition 



Tiil = TiA. 



(40) 



This condition will ensure that \E(t)\ is constant for all 
time. 

In the case of Hamiltonian systems, we would wish 
that (jl^) be a canonical transform. Since the determi- 
nant of the Jacobian of any canonical transform (e.g. E 
here) must be +1 everywhere, we are led to specify that 
\E\ is constant. [It is also required that E be symplec- 
tic, but it is well known that ([l3]) stays symplectic if 
E(to) is symplectic fT^] .] Therefore, Hamiltonian sys- 
tems demand that \E(t)\ be constant. In the case of 
dissipative systems, we still must face the boundedness 
requirements (|l4j) on Lyapunov transforms. If (^) is 
true, then \E(tf)\ — |-E(*o)|, and since under reasonable 
assumptions ( |38| ) is continuous and bounded, then |-E(t)| 
is bounded away from infinity. Finally, there is the prac- 
tical point that it is hard to program "boundedness" in 
an algorithm, but relatively easy to program constancy. 
We will make this latter choice, and enforce (|40|). 

It is now possible to explicitly determine all of the par- 
tial derivatives of the stability exponents about an equilib- 
rium point. Beginning with (fie)), we remember that each 
Qi is a function only of yi to find 



d_ 

dyi 



Trfi 



dyi 
dA aa 

dXf3 



0A o 



dyi 
Epi. 



(41) 



It is remarkable that one constraint ( ff0|) combined with 
the separation condition determines all of the local sta- 
bility exponents individually. We also note that the 
above form preserves the usual constant-coefficient case: 
dVli/dyi = if the system really is a constant coefficient 
linear problem, dAij/dxk — 0. Similarly, at the next 
order we have 



d 2 fli d 2 A 

( V ex. 



dyi 



dyf 
d 2 A aa 
dxpdx. 



-Ep; t E 1 i + 



dA aa dE, 



0i 



dxp dy l 



(42) 



This immediately generalizes to any order. Knowing the 
partial derivatives of il, it is then generally possible to 
solve (|U), (^5|), and subsequent systems of linear equa- 
tions for the coefficients of the coordinate frame trans- 
formation. Actually, we have used (24) and ( p5| ) for nu- 
merical checks. Appendix A presents an alternate form 
of these conditions that is considerably more efficient nu- 
merically. 

The one consistent case where a solution is not possible 
is where the zero order stability exponents occur as pos- 
itive / negative pairs. Such equlibrium points (including 
all Hamiltonian equilibria) produce a singular matrix at 
the third order in ( |25| ) . The question of the existence of 
separation transformations near centers and saddles will 
be investigated by other means later in this paper. 

Other choices for the stability exponents are possible, 
and we have investigated some of them. It is possible to 

2 

extremalize dQi/dyi subject to the separation condi- 
tions (p4|), in an attempt to pick maximal stability ex- 
ponents. Unfortunately, the result is always zero for the 
stability exponent partial derivatives, which certainly is 
an extremal answer. Norm-like quantities can also be ex- 
tremalized subject to the separation conditions, for ex- 
ample J^ijk (dEij /dyk) 2 is one such quantity. We also 
note that in structural mechanics there is a theory of 
"nonlinear normal modes" , Vakakis et. al ||, which min- 
imizes the curvature of the new coordinate frame. We 
have not elected that approach here. Curvature can be 
uniquely defined in the theory of structures, in general 
relativity, and in differential geometry where an under- 
lying "flat" space is implicitly assumed. But in general 
dynamical systems it is traditional to use any set of ac- 
ceptable coordinates for a system, and there may be no 
good answer to the question of which set of original co- 
ordinates is the "flat" one. 

The author believes that the choice of constant deter- 
minant of the Jacobian matrix E is compelling. It is the 
only choice for Hamiltonian systems which makes it pos- 
sible to have a canonical transformation X — * Y. For 
non-Hamiltonian systems other choices may exist which 
bound \E\, but the author is unaware of any easily speci- 
fied transformation which enforces this essential require- 
ment of the transformation. 

Since Hamiltonian systems form a good part of the 
reasons for the choice of the local stability exponents, it 
is perhaps not surprising that they assume a special form 
in this case. Partition the state vector of a 2N order 
system with Hamiltonian H as X T = {qi, p{\. Then 
direct calculation gives 



A; 



d 2 H 



A 



2 H 



Continuing, the local stability exponents are given by, in 
explicit summation notation 



(43) 



dyi 



N 2N 

EE 

a=l 13=1 



d 3 H 



dpidqidxp 



-E 



d 3 H 



Pi 



dpidqidxp 



-E, 



0i 
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(44) 



V. SEPARATION NEAR A PERIODIC ORBIT 



for Hamiltonian systems. 



IV. SEPARATION ALONG A TRAJECTORY 

The range of the cqulibrium point expansion can be 
extended by sampling the solution space along trajecto- 
ries that emanate from the equlibrium. Returning to the 
differential equation forms, we can simultaneously inte- 
grate the equations of motion (Q), the modal differential 
equations (13), and 



dt 



(45) 



To obtain a complete set of ordinary differential equa- 
tions, it is necessary to also produce a differential equa- 
tion for fij. Returning to (f4l|), and remembering that 
each Oj is a function of only one yi, we find 



4-<k 
dt 



dCli 
dyi 
OA. 



■y, 



dA n 



Ef3iE~ / 7 . 



(46) 



At this point we have a closed set of differential equations. 
In addition, 



Y — Er 1 f 



(47) 



can also be integrated to map the original coordinate 
frame X onto the new frame Y. At an equlibrium point 
all initial conditions are available: X is the equlibrium 
point state, £li = and yi = at the equlibrium point, 
and the classical decomposition furnishes the initial val- 
ues of E and f2j. Of course, a trajectory started exactly 
at the equlibrium will not evolve, but one started nearby 
will, and its initial conditions can be calculated from the 
series expansions about the equlibrium. 

Then beginning near an unstable equlibrium, we can 
integrate trajectories outward, including the stability ex- 
ponents and modal frame coordinates y. Since the modal 
transformation decouples the equations of motion and 
not the trajectories, an individual trajectory cuts across 
different coordinate lines Y, and with a dense sampling 
of trajectories the entire modal frame may be mapped 
out, starting from the vicinity of the equlibrium. 

A proof that this procedure produces an actual coor- 
dinate frame can be sketched as follows. The differential 
equations (0), @, @, @ and @ are just the differ- 
ential equations we have expanded about the equlibrium. 
When numerically integrated, they produce unique val- 
ues of Y(t) and X(£) by the standard existence theorems 
for ordinary differential equations. Then, at least locally 
X(Y) is a well defined function, since its Jacobian matrix 
<9X/9Y = E is, by construction, nonsingular. Then, as 
an immediate consequence dEij/dyk — d 2 H.i/dyjdyk is 
symmetric with respect to the indices j and fc, and Y is 
a coordinate frame. 



The construction of the standard Floquet solution (Q) 
begins with the solution to a boundary value problem 
to find periodic initial conditions. In this process, the 
variational equations (^) are integrated to help find the 
periodic orbit, and a natural by-product of this is the 
monodromy matrix $(r, 0), the state transition matrix 
at one period. Then, since the modal vector matrix E(t) 
is periodic, (fy directly shows that the eigenvectors of the 
monodromy matrix are the initial modal matrix E(0), 
while the Poincare exponents A, are related to the eigen- 
values Xi of the monodromy matrix by Aj = logAj/r. 
Then the modal matrix may be propagated for one pe- 
riod using ( ft3| ) with initial conditions E(0) known, and 
u>i — Aj taken to be constants. 

This must be modified somewhat in the current theory. 
Since the E matrix forms the basis vectors for a new 
coordinate system y, E must be periodic, and again we 
have that the matrix E(0) is the eigenvector matrix of 
3>(t, 0). But the classical Poincare exponents cannot be 
taken to be constant if the determinant of E must be 
constant. Instead, they must be interpreted through (|TJ) 
as the global stability exponents for the periodic orbit. 
This implies the constraint 



Ai = 



(48) 



relating the standard Poincare exponent to the exponents 
introduced in this work. Then, comparing known infor- 
mation to the initial conditions necessary for time propa- 
gation, we still do not know the initial values of the local 
stability exponents Qi(0). 

Taking the expression for fij ( [46| ) and integrating twice 
with respect to time along the periodic orbit gives 



fii(i-) =Qi(0)T + 



Inserting this result into (ff8|), the unknown initial condi- 
tions for the local stability exponents are 



fi i (0)=A i 



1 



dA aa 
dxp 



EpiE- 1 f 7 dt 2 



(50) 



The double integration can be easily performed with the 
time propagation algorithm by beginning the integration 
with zero initial conditions for f2j(0). There is noth- 
ing wrong analytically with (]5^) , and numerically speak- 
ing it sometimes even works. The difficulty arises when 
Poincare exponents are far from zero, leading to numer- 
ical problems in inverting an exponentially growing or 
decaying E matrix. For Hamiltonian systems this is not 
necessary, since in view of (^i|) the initial conditions are 

n,-(o) = o. 

The above suffices for isolated periodic orbits. How- 
ever, for systems with families of periodic orbit, normal 
forms for the matrix (l are required. Also, some further 



G 



considerations come into play in order to construct legit- 
imate coordinate frames Y. For two dimensional Hamil- 
tonian systems, we may begin with the usual normal form 



fi = 



1 




(51) 



for two dimensional systems. This form is only possible 
since the local stability exponents are zero for Hamilto- 
nian systems. Then, direct calculation will show that fi 
commutes with 



fi 



* 




(52) 



for any time. At the end of one period this leads to the 
extended eigenvector / eigenvalue problem 



$(t,0)E-E 



1 T 

1 



= 0. 



(53) 



As is well known, the regular eigenvector ei will be the 
state velocity vector X. The extended eigenvector will 
then be a solution of (<f> — 7)e2 = rex, combined with 
the condition ei • e2 = 0, since the starting point on an 
adjacent periodic orbit is arbitrary. 

A first integration of the differential equations of the 
previous section around the orbit will then confirm that 
the matrix E closes on itself at the end of one period, as 
it must if Y is a coordinate frame. This integration will 
also furnish the value of y\ at one period. Since ei(t) is 
the state velocity vector, the new coordinate y\ will be 
measured along the orbit itself, and herein lies a prob- 
lem. Since there are adjacent periodic orbits, and since 
the coordinate y\ must have a branch cut, it is neces- 
sary to normalize y\(r) to some constant value, in order 
that yi coordinate space not appear and disappear at 
the branch cut as we move from one periodic orbit to 
an adjacent orbit. If d\ is the multiplicative scale fac- 
tor for renormalizing ei, y\(j) is the current maximum 
value, and say 2ir is the desired maximum value, then 
d\ = 27r/j/i(r). Then, to symplectically normalize the 
E matrix, [fL2[ , the multiplicative factor di for e.^ must 
be chosen as d\di = (E T ZE) 12 , where Z is the usual 
symplectic matrix. This done, the renormalized E ma- 
trix will be symplectic, and the transformation between 
X and Y will be canonical. The branch curve for the 
coordinate yi will be the locus of starting points for the 
family of periodic orbits, and moving from periodic orbit 
to periodic orbit, the j/2 coordinate obeys 



dY 
dX 



E~ 



(54) 



which is found by differentiating (|T|) . This may be inte- 
grated to track the yi coordinate. 

The non-canonical case is somewhat more difficult, 
since the diagonal fi^ are not necessarily zero. Again 
assuming a second order system, the diagonal entries are 
given by ([ll]). For a second order system, we write 



n = 



fill S2i2 

o fi 



22 



$i= /„ n * = 1 n o 11 £ >■ ^ 



Then, we must ensure that fi and fi commute in order 
that (|lTJ) is valid. Direct calculation of fifi — fifi shows 
that this occurs if 



(fin — ^22) fii2 — ( fin — ^22 ) fil2 



(56) 



The existence of such normal forms remains a current 
research topic. 



VI. NUMERICAL EXPERIMENTS 

One system we have studied numerically is Van der 
Pol's equation 



Xl = X2, x 2 



-e{x 2 - l)x 



(57) 



for e = l. The parameters of the equlibrium point decom- 
position are listed in Table I. The origin is an unstable 
spiral point, and there is the usual stable limit cycle. This 
system makes it possible to do numerical experiments on 
both the equlibrium and periodic orbit. 



1 - 
- 

-1 - 

-2 - 




FIG. 1. Contours of 5Ryi and Sfr/i shown on the X plane 
for Van der Pol's equation. 

By integrating numerous trajectories outward from the 
vicinity of the equlibrium, and keeping track of the yi co- 
ordinate values crossed, Fig. |l| is produced. It shows the 
y frame coordinate grid with a contour separation of 0.1, 
and since the two coordinates y\ and yi are always com- 
plex conjugate for real valued x, we have shown contours 
of Sftyi and 3yi- The limit cycle is shown on the out- 
side for comparison. [Some breaks in the contours are 
numerical artefacts.] The origin of the y frame is at the 
equlibrium, and it can be seen that the coordinate frame 
is collapsing as the limit cycle is approached. This is, 
however, done maintaining a constant determinant for 
the E matrix. Apparently any coordinate frame tied to 
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the limit cycle will not match smoothly with the modal 
coordinate frame tied to the equlibrium point at the ori- 
gin. 



3 r- 

2 - 

1 - 

- 

-1 - 

-2 - 

-3 - 




-4 -3 



-1 1 
31 y, 



result, with a contour interval of 0.1. Near the equlibrium 
we have tli = 0.5 ± iy/3/2, the value from the constant 
coefficient system. The two local exponents remain com- 
plex conjugate across the X plane. The plot terminates 
on the limit cycle, since no trajectories from the equlib- 
rium can cross this line. 



FIG. 2. Van der Pol trajectories in the Kyi, Sfr/i plane. 




Another interesting way to show the trajectories is to 
plot them in the y frame itself, Fig. |. Since the plot 
again shows the real part 5iyi versus the imaginary part 
Syyi, the trajectories should not really cross each other, 
even in the y frame. To the extent that they do, we 
have pushed the numerical calculation too far in an at- 
tempt to discover what the limit cycle maps to in the 
new coordinate frame. The problem occurs because the 
individual vectors of the E matrix grow at nearly the cor- 
rect rate that the yi become constant as the limit cycle 
is approached. 




FIG. 4. Contours of 5Ryi and Sfr/i for the damped pendu- 
lum. 

The equlibrium point decomposition fails when stabil- 
ity exponents occur as positive / negative pairs, includ- 
ing centers, saddles, and in particular all Hamiltonian 
equlibrium points. We have used the damped pendulum 



x 2 , x 2 



-ex 



SIM 







(58) 



FIG. 3. Contour plot of fi on the X plane. 



as a way to explore this problem. Figure ^ shows the 
y frame [with a contour interval of 0.2] for the pendu- 
lum with a damping factor of c = 0.25. The y frame is 
relatively flat near the equlibrium point, but becomes 
severely distorted as the separatrices are approached. 
As the damping factor is decreased towards zero, this 
twisting of the coordinate frame becomes more and more 
severe near the origin, and the zone of validity of the 
linearization about the equlibrium becomes smaller and 
smaller. Then, as c — > 0, the twist of the coordi- 
nate frame tied to the equlibrium point becomes infinite. 
However, this is just what is observed numerically with 
( p5| ) [but not fl24|)] as the damping vanishes. Apparently 
this is structural, since it occurs for any choice of local 
stability exponents that the author has investigated. So, 
it appears to be impossible to use the equlibrium point 
decomposition to define a global coordinate frame when 
the equlibrium is a center or a saddle. 



A similar technique can be used to examine the behav- 
ior of the local stability exponents f2j(X. Sampling tra- 
jectories emulating from the equlibrium point, we have 



kept track of where certain values of 
then combined these into contours. 



a, 



are crossed, and 



Figure a shows the 
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FIG. 5. Three y coordinate frames based on the periodic 
orbits of the undamped pendulum. 

This does not mean that coordinate decompositions do 
not exist for the undamped pendulum. Figure || shows 
three y frames based on periodic orbits which are canoni- 
cal, and between them cover the entire phase plane. The 
yi directions are along the periodic orbits. The construc- 
tion was made with the maximum value of y\ as 2ir, since 
clearly y\ topologically looks like an angle near the stable 
equlibrium. The stable equlibrium is obviously a singu- 
lar point for this coordinate frame, and furthermore the 
saddle point at x — ±7r is also a singular point for all 
three y coordinate frames. 

We have sought numerical confirmation that using the 
methods of section IV actually produces a decoupling 
coordinate transformation. In particular, we have at- 
tempted to calculate the matrix dgi/dyj, the gradient 
vector d \E\ jdyi and dEij/dyk- The first should be di- 
agonal, with non-zero entries dgi/dyi — f2j, the second 
should be identically zero, while the last quantity should 
be symmetric in j and k if Y is a true coordinate frame. 
To calculate a numerical partial derivative, we have used 
Lagrangian five point numerical derivatives. Each nu- 
merical partial requires integrating a sheaf of five trajec- 
tories from the vicinity of the equlibrium to straddle the 
current point. Since it is common to use equal spacing in 
these points, we have solved a boundary value problem 
to find initial x values that produce equal spacing in the 
Hi coordinate. As we have chosen to make these displace- 
ments real, this means that it is necessary to integrate 
the physical system off of the real x axis. Explicitly, we 
have plotted the quantities 



a = \\E(t)\-\E(t )\\, 



(59) 



the error in propagating the determinant of E, the max- 
imum deviation of its gradient from zero, 



b = Max, 



d\E\ 



dyi 



(60) 



the maximum error in the local diagonalization of the 
equations of motion 



c = Max; 



(61) 



and the maximum violation of the symmetry condition 
for a coordinate frame 



d = Max 



]k 



dEi. 



dE, 



ik 



dyk dyj 



(62) 



as functions of time. Since these integrations extend from 
a point very close to the equlibrium well into the non- 
linear regieme, and since solution of a boundary value 
problem is required to obtain the points necessary for 
the numerical partial derivatives, the author feels that 
some numerical error growth is inevitable. 




FIG. 6. Numerical checks of the decoupling coordinate 
frame for the damped pendulum. 

Figure || shows the result of one such calculation for 
the damped pendulum with c = 0.25, a starting posi- 
tion within 0.001 radian of the origin, and whose ending 
position is well over 4 radians. After some initial error 
growth, all of these tests indicate that the time propaga- 
tion method is indeed producing a true coordinate frame 
Y which separates the equations of motion. Similar re- 
sults are shown in Fig. [7] for the Van der Pol system. 
Again the starting point is within 0.0001 units of the 
equlibrium, and the final point is quite close to the limit 
cycle. In this case error growth in the calculation of these 
four quantities is more pronounced than in the case of 
the pendulum. But the results confirm that the yi do 
decouple the equations of motion, and that the integra- 
tion is conserving the determinant \E\. Confirmation of 
the symmetry condition is the worst numerical verifica- 
tion, since this quantity is the difference of two numerical 
partial derivatives. There is good theoretical reason to 
believe that each of these conditions arc met. 
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2 4 6 8 10 12 14 16 
time 

FIG. 7. Numerical checks of the decoupling coordinate 
frame for Van der Pol's equation. 
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In this paper we have shown that the problem of local 
stability exponents which integrate along a trajectory to 
give global stability exponents is fully equivalent to the 
problem of dynamical separation of variables. The choice 
of constant determinant for the modal matrix E seems 
compelling for dissipative systems, and is the only choice 
permitted for Hamiltonian systems. The local stability 
exponents introduced by this choice are invariant to any 
linear change of variables. This invariance is less than the 
total invariance familiar from constant coefficient linear 
systems and time-periodic systems, but is still broader 
than the degree of invariance permitted by standard Lya- 
punov exponents. We have presented methods which lead 
to decoupling transformations in the vicinity of equlib- 
rium points, periodic orbits, and along general trajecto- 
ries eminating from equlibrium points. Such decoupling 
transformations may not exist in the vicinity of center 
and saddle point equlibria when the stability exponents 
of the equlibrium exist as positive/negative pairs. How- 
ever, in this case it appears possible to base decoupling 
coordinate transformations on the surrounding periodic 
orbits. 

Much work remains to be done. We are currently ex- 
ploring control applications, as well as continuing work 
on decoupling transformations near periodic orbits. An- 
other area of interest is the relation of this method to 
the "geometrodynamics" approach, which attempts to 
decouple the trajectories of a dynamical system, and 
not necessarily the underlying equations of motion. Our 
method is based on decoupling the equations of motion, 
but seems to become a trajectory-based decoupling for 
Hamiltonian systems. 



VIII. APPENDIX A. 

The form of the separation conditions near the equlib- 
rium suffer from the fact that the "columns" of the par- 
tial derivatives of E are coupled. We have used equations 
(^J) and for numerical checks, and have employed 
an alternate form for the solution that is more efficient 
numerically. 

Instead of beginning with the separation conditions, 
begin with the modal matrix equation of motion jl5|), 
written as 



Aj a Q>i8j a j E a 



Efia fa ■ 



(63) 



Evaluation at the equlibrium point immediately yields 
the eigenvalue / eigenvector problem for the first order 
terms. Then, taking a partial derivative, evaluating at 
the equlibrium, and simplifying there results 



A, 



dE ai 
dyk 



O 9E 3i 

ayk 
dtli 

— — dkiEji 
dyt 



— 

oyk 



" EykE a i 



(64) 



The advantage of this form is that each "column" (e.g. 
index i) is decoupled from all the others, reducing the 
order of the linear system by a factor of N, the order of 
the dynamical system. Treating the quantities dQi/dyi 
as known, we must solve N systems of N 2 linear equa- 
tions, not one system of order N 3 . Continuing to the 
third order gives 



A 



d 2 E ai 
Ja dy k dy, 
d 2 tl 



d 2 E n 



dy\ 



d 2 A oa 



dykdy 1 
E e iEskE ai 
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dA ja ( dE m dE Sk dE ai \ 

5 -5 E$i + — h a i + — Esk 

oxg \ ay k dyi dyi J 

dtl l<: dEu dfli; dEu 

<9y; c% dyi dy/ 

%r 0y r % 

This form also has the "column" separation property, 
since all of the second partials of £7 with second index 2 
can be solved for at once. 

Table I. 

Van der Pol Equlibrium Decomposition 

Order 1 

Q[ 0.50000000000000 + 0.86602540378444^ 

en 0.70710678118655 + 0.i 

e 2 i 0.35355339059327+ 0.61237243569579i 

Ct 2 0.50000000000000- 0.86602540378444i 

ei2 0.70710678118655 + 0i 

e 22 0.35355339059327- 0.612372435695792' 

Order 2 



fin 0. + 0.2 

em O. + O.i 

en2 0. + 0.2 

e 2 n 0. + 0.2 

e 2 i2 0. + 0.2 

fi 22 0. + 0.i 

6121 0. +0.2' 

ei22 0. + 0.2 

e 22 l 0. +0.2' 

6 2 22 0. + 0.2 



Order 3 



fim -1.00000000000000+Oi 

e im 0.0951874513135 - 0.02355278598832 

eiii2 -0.53033008588991 + 0.306186217847902 

en2i -0.53033008588991 + 0.306186217847902' 

en22 -0.53033008588991 - 0.306186217847902 

e 2 iii -0.50313367122889+0.211975073894702 

e 2 ii2 -1.0606601717798+0.2 

e 2 i2i -1.0606601717798+0.2 

e 2 i22 -1.0606601717798+0.2 

fi 22 2 - 1 .00000000000000 +0.2' 

ei2u -0.53033008588991 + 0.306186217847902' 

ei2i2 - 0.53033008588991 - 0.306186217847902' 

ei22i -0.53033008588991 - 0.306186217847902 

ei222 0.0951874513135 + 0.02355278598832 

e 22 n -1.0606601717798+0.2 

62212 -1.0606601717798+0.2 

e 2 22i -1.0606601717798+0.2 

e 22 22 - 0.503 1 33 6 7 1 2 2 8 89 - 0.21197 5 7 3 8 9 4 702' 
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